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QUANTILE ESTIMATION IN DEPENDENT SEQUENCES 



by 

P. Heidelberger P. A. W. Lewis 

IBM Research Center Naval Postgraduate School 

Yorktown Heights, NY 10598 Monterey, CA 93940 

ABSTRACT 

Standard nonparametric estimators of quantiles based on 
order statistics can be used not only when the data are i.i.d. , 
but also when the data are from a stationary, (j)-mixing process 
of continuous random variables. However, when the random vari- 
ables are highly positively correlated, sample sizes needed for 
acceptable precision in estimates of extreme quantiles are com- 
putationally unmanageable. A practical scheme is given, based 
on a maximum transformation in a two-way layout of the data, 
which reduces the sample size sufficiently to allow an experi- 
menter to obtain a point estimate of an extreme quantile. Three 
schemes are then given which lead to confidence interval esti- 
mates for the quantile. One uses a spectral analysis of the 
reduced sample. The other two, averaged group quantiles and 
nested group quantiles, are extensions of the method of batched 
means to quantile estimation. None of the schemes requires that 
the process being simulated is regenerative. 
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QUANTILE ESTIMATION IN DEPENDENT SEQUENCES 



1. INTRODUCTION 

Let {X^} be a sequence of random variables which is 

assumed to be stationary with marginal distribution 

F (x) = P{X < x} . For any p , 0 < p < 1 , the p-quantile 

X n 

of F x (x) satisfies the equation p = F x (x p ) if F x ( x ) has 

a positive density in a neighborhood of x^ . This paper is 

concerned with the estimation of quantiles x^ from samples 

X n , n=l,...,N. In particular the paper is concerned with both 

point estimates and confidence interval estimates for quantiles 

when the x n ' s are highly correlated. This is the usual case 

when the X 's are the outputs from a simulation of a stochastic 
n 

system, for example the waiting times in a queue. 

In the case of independent x n ' s > methods for quantile 

estimation are well known (e.g. Connover, 1980, pp. 71) . Thus 

let X, , denote the order statistics from the sample of size 
(n) ^ 

N , and let [xj denote the greatest integer less than or equal 
to x . Then a standard nonparametric estimator 



X p x (Lpn + lj 



is a consistent, asymptotically normal, estimator of x^ with 

bias which is 0 (N and variance [ p ( 1 — p ) ]/[ f^ (x )N] + 0(N , 

x p 

where f„ (x ) is the derivative of F w (x) at x . This vari- 
X p X p 

ance can be estimated by estimating f x (x^) < but nonparametric 

confidence intervals can also be obtained for x usinq the 

P y 

order statistics of the sample (Connover, 1980, pp. 111-116). 



2 



The difficulty with these estimation methods is that 

they require a large amount of computer storage and computing 

time to sort the sample. Also for high or low p , the order- 

statistic estimator might be biased and very non-normal. A 

solution is to use the maximum transformation (Goodman, Lewis 

and Robbins, 1971). This not only moves the problem back to one 

of estimating a median, but also allows the use of stochastic 

approximation (Robbins-Monro) methods. A very satisfactory 

method of this type using a minimal amount of memory and giving 

both point and confidence interval estimates for the unknown 

quantile x p has been given by Robinson (1975) . 

For the case of dependent X n ' s , the usual situation 

encountered in system simulation studies, quantile estimation 

is much more difficult than in the independent case. The point 

estimate x is still valid; however its variance is inflated 
P 

by a factor p(0;Xp) . Here p(0;Xp) is the initial point on 

the spectrum of the binary process {I^x^)} ' def; '- ne< ^ for each 

n to be 1 if X < x and zero otherwise. More directly it is 

n - p 

P( 0 ; x ) = lim n Var{(I (x ) +..+ I (x ))/n} . 

^ n vco F F 

The nonparametric confidence interval estimation procedures are 
no longer directly applicable. However it is possible to esti- 
mate p(0;x ) using methods of Heidelberger and Welch (1980, 

P 

1981) and to estimate the density f ( x p) by standard methods; 
used together these give an estimate of the variance of x^ 
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and large-sample confidence intervals. The main problem with 



these order-statistic estimates is , however, that the sample 
sizes required to obtain reasonable precision with the estimates 
are prohibitive when the X r ' s are highly positively correlated . 

The basic scheme considered in this paper to handle the 
sample size problem for dependent, <f>-mixing sequences is to set 
out the data in a (conceptual) v x m array 



, i X i + (k-l)m ' 

for k = 1 , . . . , v and i = 1,... m = N/v , where v is often 
but not necessarily chosen so that p V = 0.5 . Then if m(and N) 
is large enough so that X r ' s which are m apart are approxi- 
mately independent, the table can be collapsed by taking maxima 
down the columns. The maxima 



Y. 



l 



max X 
l<k <v 



for i = l,...,m are now a reduced sample whose q=p V quantile, y 

corresponds to the desired quantile, x , of the {x } process. 

p n 

The main point of this procedure is that it gives a very sub- 
stantial sample size reduction . In addition, there is generally 
slightly less overall correlation in the Y. sequence and this 
counteracts the increase in variance which occurs with the use 
of the maximum transformation and order statistic estimates of 
the quantile x^ . 



q 
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Since the Y y s (the maximum transformed sample) are 

still correlated several problems remain, notably estimating 

the variance of the point estimate y^ (the q=p V sample quantile of 

the Y-'s) and obtaining confidence interval estimates of x 
1 P 

Three confidence interval schemes are considered; averaged group 
quantiles, nested group quantiles and a scheme based on esti- 
mating the probability density function of the Y^'s at x^ , 
and the initial spectral point, p(0;y ) of the binary process 

q 

obtained from comparison of the Y^'s to y , i.e. Ip-Yq) . 

Extensive empirical sampling studies using M/G/l queues 
and stationary sequences of autocorrelated exponential random 
variables show that the above schemes provide a reliable method 
for estimating a quantile in a dependent sequence (X^i 



2. QUANTILE ESTIMATION AND THE MAXIMUM TRANSFORMATION 
1 . Quantile estimation for i.i.d. sequences . 

Let . . . ,X n ,.. . ,X N be a sample of i.i.d. random 

variables from a continuous distribution F x (x) with 
probability density function F^(x) . For 0 < p ' 1 let 

x = inf{x : F (x) = pi = F 1 (p) , 

p A A 

(p) is the inverse of F x (x) with derivative 

The quantity x is the pth quantile of F (x) 

P 



(2.1) 



where F. 



i/£ x ( y • 



Let X (1) < ,...,< x (n) £ ,. . .X (N) be the order stat- 
istics corresponding to the sample. The usual nonparametric 



point estimator of x p is the pth sample quantile 



( 2 . 2 ) 



x p X ( [Np+lJ ) ' 



0 < p < 1 



where [zj denotes the integral part of z . The following 

A 

properties of x are well known (David, 1970, 65-67) : 

P 



(2.3) 



E <V 



= x 

P 



f x ( y 

2f X (x p> 



. p . u- px + 0 

N + 2 



N 



(2.4) 



var(J ; ) = a 2 = + 0 

P P Nf x ( V 



1 



N 



and 

(2.5) 



(x - x ) N(0, No 2 ) 

P P P 



as N -*■ 00 . 

Nonparametric confidence intervals for x^ are based 
on the identity (see e.g. Connover, 1980, pp. 111-116) 



( 2 . 6 ) 



P(X (L) i X p < X (U)> = X 

j~ Ll 



N 

U 



p j d- P ) N - j , 



with L < U chosen in such a way as to make the probability as 
close as possible to the desired coverage (1-a) . 
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Another method for estimating quantiles uses the method 
of stochastic approximation (Robbins and Monro (1981)). This 
provides a method of quantile estimation which avoids sorting 
and requires only a minimal amount of storage. The asymptotic 
variance of the stochastic approximation estimate is given by 
(2.4). A very satisfactory development of this method is given 
in Robinson (1975) . Since it is not extendable to dependent 
data, it is not described. The key idea in the method, however, 
is the use of the maximum transformation, and since this is of 
use in the present context, it is described next. 

2.2. The maximum transformation 

The maximum (minimum) transformation is a computationally 
simple transformation and compaction of i.i.d. data which trans- 
forms the problem of estimating an extreme quantile of a random 
variable into that of estimating a more central quantile, e.g. 
the median. Thus assume that p > 1/2 , and let 

Y = max (X^ , . . . ,X ) . (When p < 1/2 a minimum transformation is used) 

Then 

(2.7) P { Y < x } = F (x ) = (F (x )} V = p V = q . 

— P 1 p A P 

Thus the pth quantile of is the q = p v th quantile of the 

Y variable. In particular let v = |_£n (1/2) /£n (p) J ; then x 

P 

is approximately the median of the Y variable. Details are 
given in Goodman, Lewis, and Robbins (1971) . 

There is a price paid with this transformation in infla- 
tion of the leading term of the variance (2.4). Thus let v 
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divide N , so that N/v = m . Then the X sample is reduced 



to a sample by application of the maximum transfor- 



mation to each successive non-overlapping group of X^ ' s of 

size v . Then the density of the Y.'s is f (x ) = v f v (x )p 

i y P * P 

/\. 

and the variance of the order statistic estimator, y 



v-1 



of 



in the sample Y^,...,Y is 



( 2 . 8 ) 



var (y ) 



Pd-P) 

Nf X lX p ) 



(l-p v ) 



V (1-p) p 



v-1 



var (x ) 



(1-P V ) 



p ' v(l-p)p v ‘ 1 



The right-hand multiplicative factor in this last expression is 
approximately 1.4 if v is chosen to make q = p v approximately 
0.5 . 

When using stochastic approximation, shifting to the neigh 
borhood of the median is essential for the method to be a well- 
behaved estimation procedure. For estimating several quantiles 
simultaneously, the method is applicable in a nested scheme 
which is very efficient with respect to storage and speed. 

For p < 1/2 a minimum transformation is used. Next- 
to-maximum and maximum-minimum schemes give more flexible and 
robust schemes but are not discussed in the present paper. 

2.3. Dependent Data: General Considerations and Order Statistics 

Quantile estimation in dependent data is an order of 
magnitude more difficult than in independent data. Fishman 
(1978, p. 270 ) notes that no satisfactory solution exists. 
Iglehart (1976) , Seila (1976) and Moore (1979) have given special 



8 



methods for processes {X n } with regenerative structure, i.e. 
processes for which there exist random time points at which the 
process restarts probabilistically. An example is the waiting 
time process {W } in the M/G/l queue which regenerates every 
time a customer arrives to find the queue empty, so that the 
waiting time of that customer is zero. 

It is possible to use the order statistic estimator x 

P 

given at (2.2) , where we ignore for now the problem of the 
initial transient which is usually encountered in simulations. 
Conditions for convergence and Central Limit Theorems for sample 
quantiles are well known (Sen, 1972; Babu and Singh, 1978) . 

Thus let 



(2.9) 



I (x) = 1 
n 

= 0 



X < x , 
n — 



if 

otherwise ; 



( 2 . 10 ) 



N 

I ( x ) = l I (x)/N , 
n = l 



the empirical c.d.f. at x , and 



( 2 . 11 ) 



p(0;x) = lim {N var I (x) } 

N-yco 



= l 7 v (x) , 



k=-oo 



where 

( 2 . 12 ) 



Y k (x) = COV{l n (x) , I n+k (X) 



k = 0 , + 1 , . . . 

— co ^ x < + 00 
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= P{X n i X , X n+k < x) - P{X n < x) P{X n+k < x} . 



The notation p(0;x ) comes from the fact that this 

P 

quantity is the initial point ( f =0 ) of the spectrum of the 
process {^(x^)} ; 



(2.13) 



p ( f ; x) = l cos(2TTfk) Y k (x) , 

k=-oo 



- j < f < 1/2 



Now the Central Limit Theorem for x^ = x ( | np +iJ ) 



(Sen, 1972) states that if f^(x) i s cont inuous , 

bounded and non-zero in some small neighborhood of x , and if 

I 

f^(x) is bounded in this neighborhood, then 



(2.14) 



N 1//2 (x - x ) 

P P 

(p (0 ; x ) ) 1//2 /f (x ) 
P A P 



N ( 0 , 1) 



if the process 4>-mixing and 

(2.15) o 2 = p(0;x )/{f (x ) } 2 

P P P 

is finite. For details on <J> -mixing see Sen (1972) and 
Billingsley (1968); some discussion is given in Section 2.4. 

Regenerative processes such as the M/G/l queuing 
system waiting time are ^“mixing. 

The problem with using x^ as an estimator in positively 
correlated sequences is that the sample sizes required for 
adequate precision are prohibitively large. Both sorting times 
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and memory times are then unrealistic. A measure of the infla- 
tion of sample size over the independence case is the ratio of 



identical marginal distributions. This has been investigated 
by Blomqvist (1967) for the M/G/l queue. For extreme quantiles 
of the M/M/1 queue with traffic intensity p = 0.9 this ratio 
is 400 . Greater ratios are possible, depending on the traffic 
intensity and the skewness of the service time distribution. 
Specifically, the M/M/1 queue with p = 0.9 requires a sample 
size of roughly 500,000 customers to estimate the 0.99 quantile 
of the waiting time distribution to within plus or minus 10% 
accuracy. This is derived from Table 8 of Blomqvist (1967); more 
precisely this is the sample size required for a 90% confidence 
interval for x gg to have a relative half-width of 0.10 . For 
the 0.999 quantile the required sample size is approximately 
2,300,000. Clearly storing and sorting the entire sequence is 
impractical in such cases. Actually to produce an order-statistic 
point estimate of requires storing only the largest (l-p)N 

values of the sequence. However, this ordering must be dynam- 
ically maintained as the sequence is generated, a computationally 
expensive operation. Additional storage is required to estimate 
the variance term p(0;x^) : the positions in the sequence at 

which the (l-p)N largest values occur must also be saved in 
order to construct the binary sequence (I n (Xp)} . Furthermore, 
the Discrete Fourier Transform of the extremely long sequence 




to its value p(l-p) for the i.i.d. process with 



{I (x ) , n = 1 

n p 



r • • • r 



N} must be taken to accomplish the variance 



estimation . 
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Some of the storage and sorting problems could be re- 
lieved by using stochastic approximation . However properties 
of the stochastic approximation, particularly its asymptotic 
variance, are unknown for dependent data. More importantly 

application of the maximum transformation to bring the 
quantile estimation down to a problem of estimating the median 
requires independence of the s . 

We now present a method for using the maximum transfor- 
mation to achieve sample size reduction and a practical scheme 
for point and confidence interval estimates with dependent data. 

2.4. The maximum transformation in the dependent case 

The basic idea behind the use of the maximum transfor- 
mation is to combine elements of X n , n = 1,2,..., in a <J>- 
mixing process which are sufficiently far apart so as to be 
approximately independent. To define 4>-mixing let and 

GO 

M n+m be res P ectlvel y the o-fields generated by (X^; i < n) 

and {X. ; 1 > n + m} . If E.eM and E^eM , , then 

1 — l-oo 2 n+m 

{X^} is Q-mixing if for all n(-°° < n < °°) and m(>^ 1) 

(2.16) | P (E 2 I E 1 ) -P(E 2 )|< (Mm) , Mm) > 0 , 

where 1 > 4> ( 1 ) > <J>(2) >... and lim 4> (m) = 0 . 

Thus we set out the data in a vxm array, where 
m = N/v, as 
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(2.17) 



k l/*../V/ i 1,..., 



m . 



X 



k , i 



= X 



i+ (k-1) m 



We assume for the moment that elements in the m columns 
are independent; since we are assuming that v is fixed and 
the process {X } is (^-mixing this will certainly be true as 
N (and m) get very large. Now applying the maximum transform 



down the columns 


we get a 


reduced 


series Y. 

i 


, where 




(2.18) 


Y . = 

l 


max X, 
l<k<v 


, i 




i = 1 


,2 , . . . , m . 


Values of 
mating the 


v required to 
median of the 


reduce the problem 
Y^'s are given in 


to that 
Table 2 


of esti- 

.1 . 






Size of 


Table 2 
v for 


. 1 

v . . 

P = 1/2 






p 0.900 


0. 950 


0. 975 


0. 980 


0.990 


0.995 


0.998 0.999 


v 6.6 


13.5 


27 . 4 


34 . 3 


69.0 


138.3 


346.2 692.8 



We emphasize here that transformations to points other 

than the median are possible, and are essential when considering 

simultaneous estimation of several quantiles or estimation of 

the median of the X 's . However in what follows transforma- 

n 

tion to the median is usually assumed. There will be an 
attendant bias because p V will never quite be 1/2 . 

The purposes of the maximum transformation are 
(i) to reduce the attendant sample size . It should be kept in mind 
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that for higher quantiles larger sample sizes are needed for 



given precision. The result is that for fixed precision the 
series is roughly the same length for all p > 0 . In 
particular for p = 0.999 the reduction of the sample size is 
by a factor of 693 which is sufficient to reduce a series of 
completely unmanageable length to one which can clearly be 
accommodated in a digital computer (see Section 5). 

(ii) To reduce the problem to estimation of a median rather 

than an extre me quantile . Since stochastic approximation is 
not used with the series, this is not essential for point 

estimates. It is, however, helpful in obtaining confidence 
interval estimates. 

(iii) To possibly reduce correlation . It will be shown that 
the correlations in the { Y^ } series are slightly less than 
the correlations in the (X^ } series. Though a small effect, 

it is usually sufficient to offset inflation in variances due to 

use of the maximum transformation. To see this assume that the 

structure of I (x ) is Markovian with p(0;x ) ~ (l+p)/(l-p) . 

n p p 

If p = 0.99 this equals 199. If p for the Y^ process 
chopped at its median, i.e. {l n (x Q 5 )} , is reduced to only 0.98, 
then the new p(0;x ) ~ 99 . 

t' p' ~ 

Clearly if m is too small, the assumption of inde- 
pendence down the column will be violated. Some analysis of 
this effect is possible. Expanding the definition (2.18) to 
account explicitly for the parameter v , the number of X ' s , 
v apart, over which the maximum is taken, so that 

Y. (v) = max X, . , i=l,2,...,m , we have from Billingsley 

l<k^v ' 

(1968, pp. 174, 20.49) that 
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(2.19) 



I p { Y i (v) £ Xp} - p V } | <_ vcp (m) . 



Thus we want m as large as possible and, for given v , the 
difference in the probabilities goes to zero as m -> <» , by the 
definition of mixing. A slightly tighter bound is given by the 

v — 2 , 

Lemma. |p{Y i (v) < x } - p v } < 4 , (m) p £ p k 

P k=0 



for v > 2 . 



A straight forward proof of this is obtained using induction 
The point estimate used for x is 

p 



y q Y (|_mq+lJ+l) 



where q = p v . Since the Y. ' 

1 

of this estimator are given as 
proviso for some peculiarities 
process which we discuss now. 
estimator y^ are studied by s 



s are dependent the properties 
at (2.14) and (2.15) , with 

in the structure of the {Y^} 

Finite sample properties of the 
imulation in Section 5. 



2.5. The maximum-transf ormed process {Y^} • 

The maximum- transformed process {Y^} > obtained from 

the {X n } process by taking the maximum of x n ' s which are m 
apart, is not a stationary process in the same sense that the 
{X^} process is stationary. To see this consider 
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( 2 . 20 ) 



Y, = 



max 
l<k< v 



X 1 = max{X 1 ,X m+1 , 



,X (v-l)n+l* 



It will be correlated with successive ' s , and this correla- 

tion can be expected to decrease since the {X n l process is 
mixing. However, the correlation will eventually increase because 



( 2 . 21 ) 



Y = max{X ,X„ 
m m' 2m 



/X 



vm 



} 



is coupled to by the v-1 dependent pairs 

{ VW"-- fX (v-l)m' X (v-l)m+l } • Thus the fY i’ process 
is circular. However, it is not stationary in a circular sense 

because the correlation between Y^ and 



Y 2 = max ^ x 2 ,X m+2 , . . . ,X m+2 } is through the v dependent 
pairs {X, ,X 0 ) , { X m+1 , X^_j_ 2 {X, >t _ 1 ^ t ^_ l1 ,X /tt _ 1x ^ } , and i s 



1 2 ‘ 



(v-l)m+l' (v-l)m+2 



therefore different from the correlation between Y. and Y 

I m 

In a circularly stationary process these correlations would be 

the same. Towards the middle of the process, say i = m/2 , 

the lag one correlations will be the same and there will be no 

"edge effect" because, if m is large enough, the correlation 

between Y /0 and Y. , and Y and Y will be zero. 

m/2 1 m/2 m 

Without belabouring this "edge effect" and circularity 
in the process there are two main points to be made. 



(i) For finite m it will introduce bias because, for 
instance, methods using spectral techniques for estimating 
p(0;x) in Section 3 are based on assumptions of stationary. 
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Note, however, that { Y^ } is marginally stationary, so density 
estimation techniques for fy(x) are not affected. 

(ii) For m large the effect will be negligible and 
the asymptotics go through in the usual way by ignoring a fixed 
(or slowly increasing) set of ' s at the beginning and end of 
the process. The usual assumption is that the extent of this 
set is smaller than the extent of the dependence in the {Y^} 
process . 
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3. CONFIDENCE INTERVALS FOR QUANTILES IN DEPENDENT DATA 

/N 

In the previous section we gave a point estimator 
for the quantile of the marginal distribution in a stationary 
time series. In this section we consider three methods for 
generating confidence interval estimates for these quantiles. 

The first method is an extension of the spectral method of 
Heidelberger and Welch (1980; 1981). The second and third methods 
are extensions of the method of batch means (see e.g. Law and 
Carson (1979) or Fishman (1978)). These last two methods also 
give two new point estimates. 

These methods for generating confidence interval esti- 
mates do not rely upon use of the max-transform and the case 
v = 1 corresponds to working with the original time series. 
Moreover, if the max-transform is used we do not require that 
p V = 0.50, i.e. v ~ £n.5/S,n p . There are considerable compu- 
tational savings, however, from using some maximum transformation. 
This and statistical considerations will generally dictate a 
maximum transformation using a v slightly less than the median 
transform v . Another determing factor is that one generally 
wants to estimate more than one quantile. This factor will be 
considered elsewhere. 

3.1. A Spectral Method 

Let {Y^,...,Y } be the sequence of max-transformed 
variables as defined by equation (2.18). The point estimate of 

A 

x is the qth order statistic of {Y,,...,Y } , y , where 
P 1 m q 

q = p V . Let f^(x) be the stationary density function of 
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and let p(0;x) be the stationary spectral density of the in- 
dicator function process I^x) associated with {Y^ at zero 
frequency as defined by equation (2.13). Here edge effects 
because of the quasi-circularity of the process are ignored. 

For large values of m (see Sen (1972) and Babu and Singh (1978)) 
^ 2 1/2 

^ (y - x ) / (p ( 0 ; x ) /f (x ) ) has a normal distribution with 

M F F y F 

mean zero and variance one. Confidence intervals for x based 



on this Central Limit Theorem are generated by estimating both 
p(0;Xp) and f (x^) at the estimated quantile . This is the 
first method alluded to above. 

The quantity p(0;x^) is estimated using the spectral 
method of Heidelberger and Welch (1980; 1981) applied to the 
sequence { I . ( y ) , i=l,...,m} . This method uses least squares 
to fit a low order polynomial of degree d to the logarithm of 
the first K values of the averaged periodogram of {I^y )} . 

As suggested in Heidelberger and Welch (1980; 1981) we used 
K = 25 and d = 2 , although for extreme quantiles in highly 
congested queues with short run lengths, d = 3 was required 
to produce valid confidence intervals. This point is discussed 
more fully in Section 4. 

We found little or no loss in confidence interval cov- 

/N 

erage by using the estimated quantile y with the sequence 

{ I . ( y ) } rather than using the known quantile x with the 
l q p 

sequence {i^(Xp)} . However, a proof of the convergence of 

the distributional properties of the periodogram of {l^(y^)l 
appears difficult. 
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The density f (x ) of the max-trans formed variables 
Y P 

Y, .... , Y may be estimated using a standard kernel density 
1 ' ' m 1 

estimator (see Parzen (1962) and Rosenblatt (1956)). Specifically, 
let W(x) be a weighting function such that 



0 < W(x) < W < °° , 



(3.1) 



/ W(x)dx = 1 , 



lim I xW (x) | = 0 . 
x I + 00 



Let b(m) be a sequence of bandwidth constants satisfying 
lim b (m) = 0 and lim mb(m) = 00 . For any x the kernel 

m-*-°° m^-°° 

density estimate of f^(x) is defined by 



m 



(3.2) 



f (x) = — y T—j — r 
y m b (m) 



W ( (x-Y j ) /b (m) ) 



ou 

r 1/2 

Under the conditions (3.1) and if l <f>(n) < 00 , where the 

n=0 

(f> (n) 's are defined at (2.16) , then f (x) converges in prob- 

ability to fy( x ) • 

In our case we require a density estimate at the unknown 



/\ /N 



point Xp . Simple Taylor series expansions show that f (y^) 
converges in probability to f (x ) if the first k derivatives, 

y p 

(x) , of W(x) satisfy (x) | < W < 00 , if 



/ |w (k) (x) | dx < 00 and lim | xW W (x) | = 0 and if 

— OO I I ->qo 



(k) 
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b (m) — Z^m where converges in probability to a positive 

constant and 0 < c < k/2 (k+1) . If the weighting function 

W(x) does not satisfy the different ability conditions then 
more delicate arguments are required to show convergence of 

S\ A 

f (x ) (see, for example, Robinson (1975)). 

y ir 

We have experimented with two weighting functions, a 
triangular window 



M < 1 , 

| x | > 1 , 

| x | <_ tt/2 , 

| X I > tt/2 ; 

Bandwidth sequences b (m) had the form b(m) = Z m c 

m 

for various powers of 0 < c < 1/2 and random variables Z 
r m 

which ranged from constants to measures of the spread of the 
distribution of the Y^'s . We found the density estimates to 
be relatively insensitive to both the shape of the weighting 
function and the parameters of the bandwidth sequence. For 
small samples they tended to slightly overestimate f (x ) but 
converged to f (x^) for large samples. In Section 4 we report 
the results of experiments using the triangular window with 
c = 1/3 and Z^ equal to the inter-quantile range of the 

Y j' S ' 1 * e * b(m) = { Y ( L . 7 5m+l J ) “ Y ( L . 2 5m+l J ) } m ^ * 



1 - lx 



(3.3) 



W (x) = 



and a cosine window 



(3.4) 



W (x) = 



1/2 cos(x) 
0 
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3.2. A veraged and Nested Group Quantiles: 



Definition 



We now describe two methods of generating confidence 
intervals for the quantile x^ which do not require direct es- 
timation of p(0;x ) for the { Y . } sequence. These methods, 

P ^ 

which we call averaged group quantiles (aga) and nested group 

quantiles (ngq) , are extensions of the method of batch means to 

quantile estimation. Seila (1976) considered a version of agq 

without the max transform (v=l) for regenerative processes. 

Let N be the total sample size and divide {X. . . . . . X„} 

I N 

into G non-overlapping groups with mv observations in each 
group (N = Gmv) . Define 

(3.5) X £,j,k X ( £_i) m k + j+ (k-1) m for £=1,...,G, 

j=l , . . . ,m, 

k=l , . . . , v. 



The subscript £ refers to the group number, and the data in 

each group can be thought of as forming a matrix of dimension 

v by m . The first row of this matrix is formed from the 

first m observations in the group, the second row is made up 

of the second m observations in the group, etc. For any fixed 

values of £ and j the points X. . , and X . (l<k<v) 

< , J , k jL f f K» _L 

are at lag m (i.e. , separated by m-1 observations) and for 
large values of m ' , j 1' j 2' ’ ’ ' '^£ j k are a PP rox i mate ly 

independent. Define Y. . = max X. . . . 

l<k< v ; thus Y £j is the 



jth max-transformed variable in the £th group. For large 



m 



p{Y £j 1 x p } -3 = P V with l p(Y £j i x p^ - q I < <Mm) 



v as in 



Section 2. 
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FOr each 1=1 G let Y t(l) i Y U2)---i Y Mm) be 

the order statistics of (y £1 y j2 ... y ) and define 

/S. ^ 

y qj i = Y £ ( [jnq+lJ ) ' i-e - y q £ is an order statistic estimate of 

x p derived from the 2th group. We call a group quan- 

tile estimate. As in the case of batch means, for large values 

/s /\ 

of m , ^ y q i •••»yqG^ will be distributed as i.i.d. normals 
with mean x p and variance p (0 ;x p ) /mP (x ) . This suggests 

several point and interval estimates for x 



The averaged group quantile estimate is defined by 



(3.6) 



agq (x ) 



. i I 
G ii 



y 



q2 



A confidence interval for x 

P 



is formed by assuming that 



(3.7) 



/G ( agq ( x ) - x )/S(agq) 

P P 



has a students-t distribution with G-l degrees of freedom 
where 



S 



2 



(agq) 



1 

G-l 



G 



l 

i = l 




agq ( x ) ) 



2 



This confidence interval is valid 



normals with mean x and finite 

P 

The nested group quantile 



defined to be the median of y . . 

ql / 

/\ 

y < ... < y are the order 

7 q (1) - -• 7 q (G) 



lf { yql,"-'yqG } are 

variance . 

point estimate of x 



,.,y qG - i-e- if 

statistics of y 



qi / 



i.i.d. 



is 




r 
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then the ngm is defined by (assuming G is odd) 



(3.8) 



ngq(Xp) = y q <L.5G+ij) • 



An approximate 100* (1-ct) % confidence interval for from 

(2.6) is 



(3.9; 



(y q(L) , Y q(U) 



) , 



where for any L < U, a is defined as at (2.6) , by 



(3.10] 



1 V ,G, , c ,k. G-k 

1 - a = l ( k ) ( . 5) ( . 5) 

k=L 



This confidence interval is valid if {y . ...,y ) are mutually 

CJ JL f 

independent with P{y ^ x^} = 0.5 . For large samples this 

is guaranteed by the asymptotic normality of order statistics. 

As m-*-® 0 , agq(Xp) and ngq (x^) are asymptotically un- 
biased and their confidence intervals are asymptotically valid, 
since the sequence {X n } is assumed to be mixing. Notice that 
if G = 1 , then agq(x^) = ngq(x^) = y^ where y^ is defined 
by equation (2.20) , though clearly G = 1 is not reasonable if 
a confidence interval estimate is required. Furthermore, as with 
y^ , these point and interval estimates are valid if the max- 
transform is not used (v=l) in which case one is working with 
group quantiles of the original sequence { X. } . If the max- 

transform is used there is no requirement that p V = 0 . 5 . 



24 



3.3. 



Averaged and Nested Group Quantiles: Discussion 



There are a number of possible sources of error in the 
above two confidence interval procedures. First, if m is 
small then, because of dependence, P{Y,, . <^ x } will differ 

x. 3 — p 

from q = p V by a significant amount and E(y ) ^ x . Notice 

q x, p 

that with G > 1 , m for ngq(x ) and agq (x ) is smaller 

P P 

/\ 

than the m for y and therefore this source of error is 

M 

A 

more likely with ngq(x ) and agq(x ) than with y . Second, 

i p q 

/\ 

even if the y . sequence consists of i.i.d. observations, if 
m is small E (y . ) ^ x due to the bias in order statistics for 

qx. p 

small samples. This bias will be more severe at the tails of the 
distribution. This problem is again potentially more serious with 

a 

ngq ( x ) and agq (x ) than with y . Third, the random vari- 

H ir M 

ables y . . . . ,y may be correlated. All three of these fac- 

qi , qu 

tors suggest that the number of groups, G , should be as small 
as possible. The choice G = 5 is the smallest value of G 

A /A 

for which P{y < x < y } > 0.90 (for L=l, U=G=5 

J q(l) - p J q(G) 

the confidence level is 0.9375). Given that G should be small, 

•A A 

any independence tests for y^ . . . ,y ^ will have very low 

power. The sampling experiments described in Section 4 use 

G = 5 for both agq(x ) and ngq (x ) and the correlations 

J ^ p P 



measured between the y . 's have been very low even for strongly 
correlated sequences of X n ' s (for these experiments the sample 
sizes were such that m j> 200) . Any residual correlations 
between the group medians could be further reduced by deleting 
observations between groups, i.e. by discarding observations 
before starting the next group. 
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3.3. Theoretical Comparisons 



It is interesting to compare the expected confidence 

interval widths using these three methods under the assumption 

that each method is valid. Let m = N/v and m^ = N/vG = m/G 

2 

and assume (without loss of generality) that p(0;Xp)/m^f (x^) = 1 

/\ 

Let G = 5 and assume that y •, . . . ,y are i.i.d. normals 

qi , qG 

2 

with mean x and variance p(0;x )/m f (x ) = 1 . Then the 
P PLp 

expected ngq confidence interval half-width is 1.163 
(Pearson and Hartley (1966), Table 28, p. 190). Consider now 
the agq method. The t-multiplier for a 93.75% confidence inter- 
val with 4 degrees of freedom is 2.562 so, assuming an unbiased 
estimate of the standard deviation, the expected half-width is 
2. 562//5 = 1.145 . Finally, assuming the spectral method with 

K = 25 and d = 2 produces an unbiased estimate of the point 

2 1/2 — 

estimates' standard deviation, (p(0;x ) /mf (x ) / = 1//5 . 

P P 

Then the confidence intervals expected half-width is 
2.214//5 = .99, since the effective number of degrees of freedom 
for the variance estimate was shown in Table 1 of Heidelberger 
and Welch (1981) to be 7. 

Thus, given the assumptions, the spectral confidence 
intervals (with d = 2) will on the average be narrower than both 
the agq and ngq intervals and the agq intervals will be slightly 
narrower than the ngq intervals. If a cubic polynomial is used 
in the spectral method then the spectral confidence intervals can 
be shown to be somewhat wider than the agq and ngq intervals. 

These relationships are empirically verified in Section 4. 
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4. 



EMPIRICAL STUDIES 



All the quantile estimation methods given in the previous 
section are asymptotically valid under the ^-mixing assumption. 
The crux of the matter though is whether they work with sample 
sizes required for the usual precisions required in systems simu- 
lations, say 10% to 5% as measured by ratio of the confidence 
interval half-width to the point estimate. In this section we 
use empirical sampling (simulation) to study the bias and stand- 
ard deviations of the three quantile point estimators and the 
confidence interval coverages and widths of the three quantile 
interval estimators described in Section 3. The tests were con- 
ducted on stationary sequences of correlated, exponentially dis- 
tributed random variables (NEAR(l) and GNEAR(l) processes; see 
Lawrance and Lewis (1981a) and (1981b)) and on waiting time se- 
quences in heavily congested single server queues. For each 
process the 0.50, 0.90, 0.99 and 0.999 quantiles were estimated. 



4.1. Processes Simulated 



The 


NEAR ( 1 ) process 




with 


parameters a and 3 


is defined 


as follows. Let 


fE n } 


be a 


sequence of i.i.d. ex- 


ponentials 


with mean 1 and 


let 


<V 


and {K } be mutually 
n 



independent sequences of random variables for which 



P { K = 1} = 1 - 


P{K 


= (l-oOB) = 6, where 


6 = (1-3) / (1- (1-a) 3) / 


n 


n 






and P{l n = 3} 


= 1 - 


P{I = 0} = a . Set 
n 


X Q = E q ; then 



(4.1) 



X 



n 



= I X . + K E 
n n-1 n n 



n > 0 
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is a Markovian sequence of exponentially distributed random 
variables with mean one and correlation structure 



(4.2) 



p k = corr(X n ,X n+k ) 



, 0 . k k 

= (a$) = p 



k > 0 , n > 0 



The GNEAR(l) process yields a process with exponential 
marginals and alternating positive and negative correlations. 
This process is defined by Xq = Eq and 



(4.3) 



X 



= K E + I X' 
n n n n-1 



n > 0 , 



where X' = - £n(l - exp(- X , )) . For the GNEAR(l) process 

n n- 1 

the lag one correlation is 



(4.4) p 1 =(<*$) r , 

2 

where r=l-II/6 = - 0.6449 is the maximum negative cor- 
relation attainable in a bivariate exponential distribution. 
Efficient algorithms for generating the NEAR(l) and GNEAR(l) 
are given in Lawrance and Lewis (1981b). 

The second class of processes we considered were 
(stationary) waiting time sequences in M/G/l queues, where 
the service time distribution has a hyperexponential distribu- 
tion. Specifically , let denote the waiting time of the 

nth customer and let {A , n > 1} and {S , n > 0} be the 

n — n — 

i.i.d. sequences of interarrival times and service times re- 
spectively. We assume that 
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(4.5) p { A n - = 1 “ ex P ( " Ax) ' x > 0 , A > L 

(4.6) P{S n < x} = n x {l - exp ( - v^x) } + n 2 {l - exp (- y 2 x) } , 



where y^ and are greater than or equal to zero, and 

II 1 + n 2 = 1 . Let y = (Il^/y^ + n 2 /y 2 ) 1 and let p = A/y < 1 
be the traffic intensity. Let Wq have the stationary waiting 
time distribution (see Kleinrock (1975) , p. 205; this distribu- 
tion is a probabilistic mixture of an atom at 0 and two ex- 
ponentials); then the waiting time sequence {W n l is defined 
by 



(4.7) 



W , , = (W + S - A , ) 
n + 1 n n n+1 



n > 0 



where x = max (x,0) . 

We considered five processes in testing the estimators. 
They are : 

a = 3 = 0.95; 
a = 3 = 0.995; 
a = 3 = 0.995; 



(i) NEAR ( 1 ) 
(ii) NEAR ( 1 ) 

( iii ) GNEAR(l) 
(iv) M/M/1 



(v) M/G/l 



A = 9, y x = 10, n x = 1, P = 0.90; 

A = 9, y, = 2, y 0 = 18, n. = 0.10, p = 0.90; 



(the squared coefficient of variation of the service 
time distribution is 4.556). 
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The correlation structures of these processes range from moderate 
positive correlation (NEAR(l) , a = 3 = 0.95, Pj^ ~ .9 ) to extreme 
positive correlation (NEAR ( 1 ) , a = 3 = 0.995,p R ~ .99 k ; M/M/1 
and M/G/l) to extreme negative correlation (GNEAR(l)). These 
processes were all simulated using the LL RAN DOM II random number 
generating package (see Lewis and Uribe (1981) ) . 



4.2. Simulation Results 

Tables 2-9 report the results of extensive simulation 
studies of the quantile estimation procedures detailed in earlier 
Sections. The notation used in the tables is defined as follows. 
For each process, quantile level p , value v for the max- 
transform, and run length N , R i.i.d. replications were per- 
formed (R = 200 for all runs, except runs in which p = 0.999, 

and runs in which p = .99 with large values of m) . However, the 
same random number seeds were used for all processes, p, v and 

N so that different entries in the tables are not independent. 

Let y (r) , ngq (x , r) and agq (x , r) denote the realiza- 
^1 P P 

/N 

tions on the rth replication of , ngq(x^) and agq(x^) 



respectively. Let y , ngq and agq denote the averages 
over the R replications of y (r) , ngq (x ,r) and agq(x ,r) 

~1 IT P 

respectively. For example 



(4.8) 



? q = X V r ’ /R 



Let sd (y ) , sd(ngq) and sd(agq) denote the sample standard 



deviations of y^ , ngq and agq respectively. For example 
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(4.9) 



sd 2 (y a ) = I (y„(r) - y ) 2 /r(r-1 ) . 
q r=l q q 

These empirical sampling, point estimates and their re- 
spective estimated standard deviations may be used to study the 
bias in the quantile estimates derived in this paper. For 
example, from Table 2, an approximate 90% confidence interval 

/V 

for E (y ) in the NEAR ( 1 ) process with a = 6 = 0.95, p = 0.99, 

N = 69000 and v = 69 is V q ± 1.645 sd(y^) = 4.599 + 1.645 x .009 

The true value x 99 = 4.605 lies within this 90% confidence 

interval, which is (4.584, 4.614). 

(i) Table 2; point estimates of quantiles 
Table 2 illustrates the effect of the maximum transfor- 
mation on y^ and its standard deviation, sd(y^) . For each 
process, run length N and p = 0.90 and 0.99, y anc 3 sd(v ) 

q q 

were computed using v = 1 (no max transform) and v chosen 

so that p V ~ 0.50. The run lengths N were chosen so that 

m = N/V = 1000 for the NEAR ( 1 ) and GNEAR(l) processes and 

m = N/v = 2000 for the M/M/1 and M/G/l queues. These yield 

conservative values of m. The only case in which the max- 

transform introduces apparent bias is the M/M/1 queue with 

p = 0.90, v = 7, N = 14000; in this case v differs from x 

P P 

by 2.45 times the estimated standard deviation. Since this is 
the maximum deviation of 20 experiments, it is probably not 
significant . 

The column labeled "Actual sd inflation" is (with p 
and N fixed) the ratio sd(y I p V ~ 0.50)/sd(y |v=l) . This 

ratio measures the increase in standard deviation in dependent 
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sequences due to the max-transform. For example, for the NEAR ( 1) 
with a = 3 = 0.995 , p = 0.99 and N = 69000 this ratio is 
.0090/. 0087 = 1.03 . 

The column labeled "Theoretical sd inflation for i.i.d. 
samples" is the inflation in standard deviation due to the max- 
transform for a sequence if i.i.d. random variables; see 
equation (2.8). Recall that the max-transform changes the cor- 
relation structure of the process so that its effect on the 
standard deviation in dependent sequences may differ substantially 
from its effect in independent sequences. In most cases the 
actual inflation in variance is very slight and in several cases 
a small variance reduction is achieved. Only for the M/M/1 
queue is the inflation greater than what it would be for an i.i.d. 
sequence, and even in this case it is, considering the storage 
and computational savings of the max-trans form, an acceptable 
1.35 . 

From Table 2 we conclude that, for sensibly large values 
of m , the max-transform introduces very little, if any, bias 
and that the inflation in variance is modest. 

(ii) Tables 3-8; confidence interval estimates 
Tables 3-8, in addition to reporting the estimated means 
and standard deviations of the quantile point estimates, 
compile the results of confidence interval coverage studies for 
the three quantile confidence interval procedures. For each 
process, p , v , N and replication number, confidence inter- 
vals for Xp were generated using the spectral method, nested 
group quantiles, and averaged group quantiles as described in 
Section 3. 
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For these studies, G , the number of groups was set at 
5 . For this value of G , under the assumptions of Section 3 
the ngq confidence interval has nominal coverage 0.9375, i.e. 

A /V 

P{y q(l) - X p < y q(5) J = 0.9375. For comparison purposes we also 
formed 93.75% confidence intervals using the spectral and agq 
confidence intervals. Let CI r (y g ) , CI r (ngq) and CI r (agq) 
denote the confidence interval on the rth replication using 
the spectral , ngq and agq methods respectively and let 

A 

| Cl r ( y ) | , |ci r (ngq) | and |CI r (agq) | denote the widths of 

these intervals. For each confidence interval procedure 

Tables 3-8 report the estimated average confidence interval 

relative half-width (labeled hw/x ) . For example, for the 

P 

spectral confidence intervals 

R 

(4.10) hw/x = ( 1/R) l |CI (y )|/2x . 

F i=l q p 

These tables also report the fraction of these confidence 

intervals which actually contain x . This fraction is called 

P 

an (estimated) 93.75% coverage and it should be close to 
0.9375 if valid confidence intervals are being formed. A one 
sided confidence interval for a coverage, based on the normal 
approximation to the binomial distribution, can be used to test 
the hypothesis that the actual coverage is less than 0.9375 . 

For R = 200 (Tables 3, 4, 5, 6 and part of 8) estimated coverages 
less than or equal to 0.91 are significantly low at the 0.90 
level while for R = 100 (Table 7 and part of 8) estimated 
coverages less than or equal to 0.89 are significantly low at 
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the 0.90 level. However, in practice, the importance of cover- 
age in judging the quality of a confidence interval procedure 
depends very much on the accuracy (as measured by relative half- 
width) of the confidence intervals. Thus low coverage is not a 
serious drawback if the procedure generates confidence intervals 
which are very wide (very inaccurate) and therefore provide very 

limited information about the quantity being estimated. In fact, 

♦ 

the most useful information such a confidence interval often 
provides is that the run length is too short. As the relative 
half-width decreases, the importance of coverage for the va- 
lidity of a procedure increases. We note that it would also be 
possible to plot the coverage functions for the spectral and 
agq confidence intervals (Schruben (1980)); however, reporting 
the individual 93.75% coverages is representative of the cover- 
age function and more compact. 

For the spectral method and the M/M/l and M/G/l queues 
the coverages using both degrees d = 2 and d = 3 are given 
whereas for the NEAR(l) and GNEAR(l) processes only d = 2 is 
given. For example, from Table 4, the coverages for in the 

M/M/l queue with p = 0.90, v = 1 and N = 14000 for the 
spectral method using degrees d = 2 and 3 , nested group 
quantiles and averaged group quantiles are .905, .910, .880 

and .870 respectively. The average relative half-widths are 
.549, .570, .429 and .403 respectively. Thus to increase the 

precision from its current value of about 50% to a desired value 
of 10%, we estimate that the sample size would have to be in- 
creased by a factor of 25 to a total of 350,000 observations. 
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It would clearly be impractical to store and sort all of these 
observations. Thus to increase the precision to an acceptable 
level requires the max-transform to reduce the sample to a 
manageable size. 

Notice that when the coverages of all methods are valid 
(for example in Table 5 for all run lengths of the NEAR ( 1 ) and 
GNEAR(l) processes and large run lengths of the queueing processes) 
that the relative widths of the confidence intervals are generally 
as predicted in Section 3, namely that 

hw(spectral, d=2) hw(agq) £ hw(ngq) hw(spectral, d=3) . 

Tables 3 and 4 list results of experiments for estimat- 
ing the median and 0.90 quantile respectively without using the 
max-transform (v = 1). The run lengths used were small enough 
to accommodate storing and sorting the full sequences. 

In Table 3 all coverages are either acceptable (i.e., 
not significantly less than 0.9375) or nearly acceptable (i.e., 
wide confidence intervals with coverage close to, but signifi- 
cantly less than 0.9375) with the exception the GNEAR(l) process 

using the spectral method. The reason for the low coverage in 

/\ 

this case is that the sequence of indicator functions {l^(x <-) } 
is very nearly a deterministic sequence of alternating zeros 
and ones due to the strong alternating negative and positive 
correlation structure of the GNEAR(l) process. This phenomenon 
appears to be peculiar to the median (the coverage for the 0.90 
quantile using the spectral method in Table 4 is acceptable) but 
should be kept in mind when dealing with processes having this 
type of correlation structure. 
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The coverages for the 0.90 quantiles in Table 4 are all 
acceptable or nearly acceptable. However, agq and ngq exhibit 
substantial bias in the M/G/l queue. In both Tables 3 and 4 
notice the very large mean relative half-widths in the M/M/1 
and M/G/l queues, again implying the need for the max-transform 
to deal with the much larger sample sizes required for estimates 
of acceptable precision. 

Tables 5, 6 and 7 compile results of experiments for 
estimating the 0.90, 0.99 and 0.999 quantiles, respectively, 
using the max-transform with v chosen so that p ~ 0.50 . 

The ngq and agq confidence interval estimate coverages 
are generally acceptable throughout these tables. However, as 
in the case of v = 1 , ngq and agq show that there is sub- 
stantial bias in the ngq and agq point estimates in the 
M/G/l queue for small N. For N = 224,000 this bias has dis- 
appeared even though the precision (about 20%) is still relatively 
high. The estimates x^ show little, if any, bias. The cover- 
ages of the confidence interval estimates in the M/G/l queue 
are high despite the bias due to the extreme width of the confi- 
dence intervals at the small run lengths. 

For the spectral method the coverages for the NEAR(l) 
and GNEAR(l) processes are acceptable. However, the small sample 
coverages for the queues, particularly the M/G/l queue, are 
disappointing. Even using a cubic polynomial to fit the log of 
the periodogram does not yield acceptable coverages for the 
M/G/l queue when N is small. The low coverage is explained 
by the fact that the spectrum of the max-transformed variables 
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mv , 



{ Y k } is still very steep near f = 0 (i.e., the Y^s are 

highly correlated) even for large values of v . Since N = 
if v is large then m is small and the least squares fits to 
the log of the periodogram is done over a relatively large fre- 
quency range. If log (p(f ; x ) ) is very steep then it may not 
be well approximated by a low order polynomial over a wide fre- 
quency range. 

Figure 1 illustrates this problem in the case of the 
NEAR ( 1 ) process with a = 3 = 0.995 . The log of the spectrum 
of the max-transformed variables (v = 693) is seen to be still 
quite steep near f = 0 . Compare this to the almost flat log 
of p (f) / the spectrum of the process of batch means X_(k) , 

where 



(4.11) 



V k: 



1 

B 



kB 

I 

j= (k-1 ) B+l 



X 



These batch means would be used to place a confidence interval 
on E(Xj) (see Heidelberger and Welch (1981)). 

The problem can be avoided by making v smaller, and 
thus m larger. The least squares fit is then done over a 
narrower frequency range and the coverage improves. This effect 
is seen by comparing the coverages for x in the M/G/l 

queue with N = 14000 and v = 1 (Table 4) to those using v = 7 
and N = 14000 (Table 5) . (Again note that the precision is 
unacceptably high) . It is also seen in Table 8 which, for the 
M/M/1 and M/G/l queues, compares the coverages for x 
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Figure 1. Effect of the maximum transform and batching 
on the log-spectrum. Figure la gives the log 
of the spectrum of a NEAR(l) process with 
a = 3 = 0.995. Figure lb shews the log of 
the spectrum after the maximum transformation 
with v = 693. A small but significant effect 
can be seen 'at f = 0. Figure lc shows the 
log of the spectrum after batching, with batch 
sizes of v = 693. Batching reduces the initial 
point of the spectrum much more than the maximum 
transformation. 
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log spectrum of max transformed variables ^v=t>yj; 




Figure 1. Effect of the maximum transform and batching 
on the log-spectrum. Figure la gives the log 
of the spectrum of a NEAR(l) process with 
a = 3 = 0.995. Figure lb shows the log of 
the spectrum after the maximum transformation 
with v = 693. A small but significant effect 
can be seen at f = 0. Figure lc shows the 
log of the spectrum after batching, with batch 
sizes of v = 593. Batching reduces the initial 
point of the spectrum much more than the maximum 
transformation. 
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log spectrum of batch means. Batch size = 693 




Figure 1. Effect of the maximum transform and batching 
on the log-spectrum. Figure la gives the log 
of the spectrum of a NEAR(l) process with 
a = 6 = 0.995. Figure lb shews the log of 
the spectrum after the maximum transformation 
with v = 693. A small but significant effect 
can be seen at f = 0. Figure lc shows the 
log of the spectrum after batching, with batch 
sizes of v = 693. Batching reduces the initial 
point of the spectrum much more than the maximum 
transformation. 
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using v = 69 to those using v = 17 and v = 8 , and the 
coverages for *.999 using v = 693 to those using v = 173 
and v = 86 (the slight differences in N are to make m 
highly composite which increases the computational efficiency of 
the Discrete Fourier Transform used to calculate the periodogram) . 
Thus the full storage and computational savings of the max- 
transform (i.e. moving the quantile back to the median) cannot 
be realized in highly correlated sequences without sacrificing 
confidence interval coverage. Significant savings can be realized 
however. For example the max- transform with v = 173 reduces 
the 1,384,000 observations to a sample size of just 8000 for 
estimating x in the M/G/l queue. 

In Table 8 two simulations are also shown with 
N = 5,536,000 and v = 173 for both the M/M/1 and M/G/l 
queues for the extreme case of the .999 quantiles. This sample 
size is large enough to obtain almost 10% precision. To within 
the precision of the empirical sampling experiment (R = 100) , it 
should be noted that all point and confidence interval estimates 
are valid. Note that the maximum transformation reduces the 
sample size from 5,536,000 to 32,000!! 
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5 . 



SUMMARY AND CONCLUSIONS 



In this paper we have investigated point and interval 
estimates for a single quantile in a fixed length sequence of 
dependent observations. For large sample sizes it becomes im- 
practical to store and sort the entire sequence. For extreme 
quantiles, these limitations can be overcome by using the maxi- 
mum transformation which requires storing only a sequence of 
maxima. This sequence is defined by laying the data out into a 
v by m array and storing only the maximum element in each 
column. Storage requirements are thus reduced by a factor of v . 
Observations at lag m are assumed to be independent and the 
max- transform changes the problem of estimating the p quantile 
of the original sequence into one of estimating the p v quantile 
of the max- transformed sequence. 

Three confidence interval methods which can exploit the 
sample size reduction produced by the max- transforma tion were 
described and tested; the spectral method and two extensions to 
the method of batch means called nested group quantiles (ngq) 
and averaged group quantiles (agq) . The three methods produced 
valid confidence intervals if the sample sizes were large enough 
to produce 10% precision in the estimates. 

Problems which were not addressed in this paper and are 
the subject of ongoing research involve questions of multiple 
quantile estimation and the incorporation of these fixed run 
length procedures into sequential procedures for, say, simulation 
run length control. Among the issues involved here are: 
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(i) Testing the adequacy of the spacing parameter m 



so that and X^ +rn are approximately independent. 

(ii) Organizing the data so that m can be increased 
as the sample size increases. 

(iii) Organizing the data so that multiple quantiles can 
be efficiently estimated, for example, estimating a 0.50 and 
0.999 quantile for the same set of observations (this probably 
involves the use of a max-min transformation) . 

(iv) Determining the sensitivity of these schemes to 
an initial transient which is often encountered in simulation. 
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Table 2. Comparison of point estimates and standard deviations of point estimates 

using full sample order statistic and max- transformed sample order statistic 
(R = 200 replications). Here q = p v . The fact that the actual sd inflation 
is generally less than the theoretical sd inflation for i.i.d. sequences is 
because of the decrease in correlation in the series induced by the maximum 
transformation. 



Process 


P 


X 

p 


N 


V 




sd(y ) 
9 


Actual 

sd 

inflation 


Theoretical sd 
inflation for 
i.i. d. sequences” 


NEAR(l) 


0.90 


2.303 


7,000 


1 


2.304 


.0104 






a=6=0. 95 


0.90 


2.303 


7,000 


7 


2.320 


.0106 


1.02 


1.18 




0.99 


4.605 


69pOO 


i 


4.610 


.0087 








0.99 


4.605 


6^000 


69 


4.599 


.0090 


1.03 


1.20 


NEAR(l) 


0.90 


2.303 


7,000 


1 


2.344 


.0325 




| 


a= 6=0. 995 


0.90 


2.303 


7,000 


7 


2.340 


.0356 


1.10 


1.18 




0.99 


4.605 


69,000 


1 


4.540 


.0275 




| 




0.99 


4.605 


69,000 


69 


4.607 


.0299 


1.09 


1.20 


GNEAR(l) 


0.90 


2.303 


7,000 


1 


2.282 


.0132 






a=6=0. 995 


0.90 


2.303 


7,000 


7 


2.291 


.0121 


0.92 


1.18 


negative 


















correlation 


0.99 


4.605 


69P00 


1 


4.605 


.0046 




1 




0. 99 


4.605 


69,000 


69 


4.600 


.0050 


1.09 


1.20 

1 


M/M/1 


0.90 


2.197 


14,000 


1 


2.222 


.0277 




1 


p=0. 90 


0.90 


2.197 


14P00 


7 


2.289 


.0375 


1.35 


1. 18 




0.99 


4.500 


138P00 


1 


4.514 


.0368 








0.99 


4.500 


138000 


69 


4.558 


.0455 


1.24 


1.20 1 


M/G/l 


0.90 


6.311 


14,000 


1 


6.466 


.1682 






P=0 . 90 


0.90 


6.311 


14P00 


7 


6.572 


.1759 


1.05 


1.18 




0.99 


13. 130 


138P00 


1 


13.182 


.2077 








0.99 


13.130 


138,000 


69 


13.398 


.1975 


0.95 


1.20 
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Table 3. Point estimates, standard deviations of point estimates, 93.75% 
coverages, and average confidence interval half-width/x for the 
three methods. Here p = 0.50, v = 1, G = 5, and R, the^number 
of replications, is 200. 

T 

I 



1 

4 

j>cess 


X 

P 


N 


Spectral Method 


Nested Group 
Quantiles 


Averaged Group 
Quantiles 


X 

P 

sd(x ) 
P 


degree of polynomial 
3 2 




coverage 

hw/x 

P 




coverage 

hw/x 

P 


coverage 

hw/x 

P 


coverage 

hw/x 

P 


ngq 

sd(ngq) 


ngq 

sd(agq) 


IpAR(l) 

f 3= 0.95 


0.693 


7p00 


. 696 
(.003) 




.935 

(.149) 


. 700 
(.004) 


.935 

(.180) 


. 706 
(.003) 


.945 

(.167) 

j 


EAR(l) 

x 3=0.995 

111 


0.693 


7,000 


. 735 
(.012) 




.865 

(.517) 


.744 

(.015) 


.905 

(.598) 


.809 

(.013) 


.910 

(.565) 


;ear(d 

(3=0. 995 

iative 

(relation) 


0.693 


7,000 


.692 

(.001) 




. 795 
(.006) 


.693 

(.002) 


.945 

(.112) 


.689 

(.002) 


.985 

(.102) 


/M/1 

0.90 


0.588 


14P00 


.601 

(.006) 


.950 

(.452) 


.950 

(.352) 


.598 

(.007) 


.945 

(.443) 


.640 

(.008) 


.940 

(.413) 


I/G/l 
: 0. 90 


1.545 


14p00 


1.607 

(.033) 


.950 
(. 878) 


.940 

(.647) 


1.562 

(.036) 


.930 

(1.013) 


1.943 

(.055) 


.940 

(.961) 

1 

J 
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Table 4. Point estimates, standard deviations of point estimates, 93.75% 
coverages, and average confidence interval half-width/x^ for the 
three methods, no maximum transform. Here p = 0.90, v = 1, 

G = 5 and R, the number of replications, is 200. 









Spectral Method 


Nested Group 
Quantiles 


Averaged Group 
Quantiles 










degree of polynomial 
3 2 










Process 


X 

p 


N 


X 

P 

sd(x ) 
P 


coverage 
hw/ x 

P 


coverage 

hw/x 

P 


ngq 

sd(ngq) 


coverage 

hw/x 

P 


ngq 

sd(agq) 


coveragi 
hw/ x 

P 


NEAR(l) 
a=3=0. 95 


2.303 


7,000 


2.304 

(.010) 




.930 

(.129) 


2.287 

(.013) 


.925 

(.156) 


2.296 

(.010) 


.915 

(.145) 


NEAJR(l) 
a= 3=0.995 


2.303 


7)300 


2.344 

(.032) 




.935 

(.468) 


2.189 

(.038) 


.955 

(.469) 


2.264 

(.030) 


.940 

(.435) 


GNEAR(l) 
a= 3=0.995 
(negative 
correlation) 


2.303 


7,000 


2.305 

(.010) 




.925 

(.141) 


2.282 

(.013) 


.955 

(.181) 


2.267 

(.011) 


.955 

(.168) 


M/M/1 

p=0. 90 


2.197 


14,000 


2.222 

(.028) 


.910 

(.570) 


.905 

(.549) 


2.066 

(.027) 


.880 

(.429) 


2.229 

(.027) 


.870 

(.403) 


M/G/l 
p=0. 90 


6.311 


14P00 


6.466 

(.168) 


.890 

(1.069) 


.845 

(.669) 


5.266 

(.090) 


.865 

(.496) 


5.843 

(.099) 


.815 

(.466) 
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Table 5. Point estimates, standard deviations of point estimates, 93.75% 
coverages, and average confidence interval half-width/x for the 
three methods. As in Table 4, p = .90, but here v = 7; P thus p v ~ 0.50 
Also G = 5 and R, the number of replications, is 200. 



-I 

-I 

>pcess 


X 

p 


N 


Spectral Method 


Nested Group 
Quantiles 


Averaged Group 
Quantiles 


sd(y q ) 


degree of 
3 


polynomial 

2 




coverage 
hw/ x 

P 


ngq 

sd(agq) 


i 

coverage | 
hw/x 

P 

i 


coverage 

hw/x 

P 


coverage 

hw/x 

P 


ngq 

sd(ngq) 


EAR ( 1 ) 


2.303 


7,000 


2. 320 




.945 


2.310 


.940 


2. 328 


.940 


i 8=0. 95 






{.Oil) 




(.156) 


(.015) 


(.179) 


(.012) 


(.164) 






28,000 


2.304 




.940 


2. 303 


.920 


2.305 


.915 


1 






(.006) 




(.075) 


(.007) 


(.089) 


(.006) 


(.083) 






112,000 


2.300 




.930 


2.296 


.920 


2.300 


.935 








(.003) 




(.037) 


(.003) 


(.044) 


(.003) 


(.626) ! 


'IEAR(l) 


2. 303 


7,000 


2. 340 




.800 


2.253 


.940 


j 2.414 


.945 


i 3=0.995 






(. 036) 




(.342) 


(.038) 


(.519) 


! (.033) 


(.488) 


If 




28,000 


2. 311 




.935 


2.297 


.960 


2. 350 


.940 








(.017) 




(.246) 


(.020) 


(.279) 


(.017) 


(.266) | 






112,000 


2.289 




.970 


2.301 


.945 


2. 309 


.945 








(.008) 




(.127) 


(.010) 


(.142) 


(.008) 


(.131) : 


EffiAR(l) 


2.303 


7,000 


2.291 




.960 


2.316 


.9 35 


2.287 


.940 


=6-0. 995 






(.012) 




(.199) 


(.013) 


(.192) 


(.012) 


(.178) 






28,000 


2.306 




.940 


2. 312 


.940 


2.299 


.935 


egative 






(.006) 




(.078) 


(.007) 


(.667) 


(. 8(96) 


(.686) | 


orrelat ion) 


















1 






112 ,000 


2.304 




.925 


2.306 


.925 


2.304 


.935 








(.003) 




(.040) 


(.0 03) 


(.047) 


(.003) 


(.043) 


m/m/ 1 


2.197 


14,000 


2.289 


.940 


.815 


2.206 


.955 


2.331 


.930 


=0.90 






(.037) 


(.500) 


(.306) 


(.033) 


(.451) 


(.032) 


(.422) 






56,000 


2.231 


.975 


.950 


2.220 


.945 


2.316 


.935 








(. 016) 


(.325) 


(.240) 


(.020) 


(.296) 


(.020) 


(.274) 






224,000 


2.206 


.965 


.950 


2.204 


.955 


2.227 


.945 








(.008) 


(.142) 


(.116) 


(.009) 


(.238) 


(.668) 


(.127) 


M/G/l 


6.311 


14,000 


6.572 


.695 


.525 


5.657 


.865 


6.112 


.860 


=0.90 






(.176) 


(.399) 


(.219) 


(.123) 


(.500) 


(.110) 


(.222) 






56,000 


6.373 


.950 


.870 


6.184 


.940 


6.680 


.935 








(.079) 


(.471) 


(.300) 


(.079) 


(.444) 


(.085) 


(.417) 






224,000 


6.388 


.970 


.955 


6.367 


.965 


6.559 


.960 








(.042) 


(.260) 


(.209) 


(.047) 


(.243) 


(.047) 


(.227) 



Table 6. Point estimates, standard deviations of point estimates, 93.75% 
coverages, and average confidence interval half-width/x for the 
three methods. Here p = 0 . 99 , v = 69 , so that P 

q - P V £ 0.50 



is 200. 



Also G = 5 and the number of replications, 
However R = 100 if the run is marked t . 



R, 



Process 


X 

p 


N 


Spectral Method 


Nested Group 
Quantiles 


Averaged Grouj 
Quantiles 


sd(y q ) 


degree of 
3 


polynomial 

2 




coverage 

hw/x 

P 




\ 

cove r< 
hw/ x 

I 


coverage 

hw/x 

P 


coverage 

hw/x 

P 


ngq 

sd(ngq) 


ngq 

sd(agq) 


NEAR(l) 


4.605 


69,000 


4.599 




.960 


4.606 


.920 


4.627 


. 91! 


a=6=0. 95 






(. 009 ) 




(.063) 


(.011) 


(. 073 ) 


(.009) 


(.067 






138,000 


4.608 




.945 


4.604 


.950 


4.608 


. 95C 








(.006) 




(.043) 


(.007) 


(. 054 ) 


(.006) 


(.05C 






276,000 


4.604 




.930 


4.603 


.945 


4.608 


. 95C 








(.004) 




(.031) 


(.005) 


(. 037 ) 


(.005) 


(.934 


NEAR(l) 


4.605 


69,000 


4.607 




.890 


4.606 


.965 


4.692 


.965 








(. 050) 




(.180) 


(.032) 


(. 248) 


(.529) 


( . 225 






138,000 


4.653 




.930 


4.653 


.905 


4.685 


. 925 








(.023) 




(.140) 


(.027) 


(.173) 


(. 024) 


(. 256 






276,000 


4.620 




. 940 


4.613 


.945 


4.629 


.935 








(.014) 




(.099) 


(.016) 


(.111) 


(.014) 


(.10k 


GNEAR(l) 


4.605 


69,000 


4.600 




.955 


4.596 


.950 


4.594 


. 94C 


ot=B=0. 995 






(.005) 




(.035) 


(.006) 


(.040) 


(.995) 


(.93? 






138,000 


4.605 




.920 


4.605 


.925 


4.601 


. 9 1C 


(negative 






(.004) 




(.023) 


(. 005 ) 


(.028) 


(.004) 


(.926 


correlation) 




276,000 


4.608 




.950 


4.604 


.930 


4.605 


.935 


1 

i 






(.003) 




(.017) 


(.003) 


(.020) 


(.995) 


(.015 


M/M/1 


4.500 


138,000 


4.558 


.865 


.730 


4.351 


.920 


4.538 


.905 








(.045) 


(. 262) 


(.153) 


(. 035 ) 


(.267) 


(.955) 


(.259 






276,000 


4.513 


.955 


.895 


4.472 


.930 


4.598 


. 93C 








(. 027) 


(. 236) 


(.148) 


(. 032 ) 


(.222) 


(.925) 


(.207 






552,000 


4.525 


.955 


.945 


4.511 


.945 


4.608 


.935 








(.020) 


(.190) 


(.132) 


(. 024 ) 


(.178) 


(.925) 


(.166 


t 




2,208,000 


4.481 


.980 


.960 


4.495 


.960 


4.510 


.950 








(. 012 ) 


(.084) 


(.067) 


(. 015 ) 


(.076) 


(.013) 


(.070 


M/G/l 


13.130 


138,000 


13.398 


.675 


.490 


11. 713 


.820 


12.174 


.815 








(.198) 


(. 243 ) 


(.129) 


(. 125 ) 


(.255) 


(.119) 


(.279 






276,000 


13.359 


. 845 


.635 


12. 377 


.920 


12.934 


.900 








(. 164 ) 


(.270) 


(.155) 


(. 107) 


(.287) 


(.106) 


(.257 






552,000 


13.271 


.935 


.815 


12.927 


.910 


13.480 


.915 








(. 106) 


(. 276) 


(.164) 


(. 103 ) 


(.265) 


(.101) 


(.247 


+ 




2,208,000 


13.102 


. 960 


.930 


13.010 


.920 


13.241 


.900 








(.062) 


(.163) 


(.123) 


(. 075 ) 


(.133) 


(. 066) 


(. 225 



Table 7. 



Point estimates, standard deviations of point estimates, 93.75% 
coverages, and average confidence interval half-width/x for the 
three methods. Here p = 0.999, v = 693, so that p 

9 = P V ^ 0.50 . Also G = 5 and the number of replications, R, 
is 100. 
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Table 8. M/M/1 queue. Point estimates, standard deviations of point 

estimates, 93.75% coverages, and average confidence interval 
half-widths/x comparing q = p v ~ 0.50, q = p v ~ 0.84 and 
q = pV ~ 0.92 P for different sample sizes. Here p = 0.99 and 
p = 0.999 and G = 5. The number of replications is R = 200 
for runs marked * , and R = 100 for runs marked t . 











Spectral Method 


Nested Group 
Quantiles 


Averaged Group 
Quanti les 












degree of 
3 


polynomial 

2 














P 




y q 


coverage 


coverage 


ngq 


coverage 


ngq 


coverag 


Process 




X 

P 


N 


sd(y q ) 


hw/x 

P 


hw/x 

P 


sd(ngq) 


hw/x 

P 


sd(agq) 


hw/x 

P 


M/M/ 1 

p=0. 90 


& 


0.99 

4.500 


138,000 

69 


4.558 

(. 045 ) 


.865 

(. 262 ) 


. 730 
(. 153 ) 


4.351 

(. 035 ) 


.920 
(. 267 ) 


4.538 

(. 035 ) 


.905 

(. 250 ) 




* 


0. 99 
4.500 


136,000 

17 


4.546 
(. 040 ) 


.925 

(. 375 ) 


. 880 
(. 246 ) 


4.258 

(. 029 ) 


.895 
(. 240 ) 


4.409 

(. 029 ) 


.890 

(.222) 




* 


0.99 

4.500 


128,000 

8 


4.596 

(. 042 ) 


.900 

(. 384 ) 


.870 

(. 306 ) 


4.317 

(. 035 ) 


.895 

(. 260 ) 


4.470 

(. 037 ) 


.875 

(. 243 ) 




t 


0.999 

6.800 


1,386,000 

693 


6.835 

(. 058 ) 


.900 

(. 181 ) 


. 730 
(. 108 ) 


6.648 

(. 046 ) 


.930 

(. 172 ) 


6.790 

(. 046 ) 


.940 

(. 160 ) 




t 


0.999 
6. 800 


1,384,000 

173 


6.855 
(. 059 ) 


.910 

(. 284 ) 


.870 

(. 178 ) 


6.578 

(. 055 ) 


.930 

(. 178 ) 


6.732 

(. 050 ) 


.910 

(. 165 ) 




t 


0.999 

6.800 


1,376,000 

86 


6.856 
(. 056 ) 


.940 

(. 232 ) 


.890 

(. 181 ) 


6.657 

(. 054 ) 


.900 

(. 156 ) 


6.754 
(. 049 ) 


.880 

(. 144 ) 




t 


0.999 

6.800 


5,536,000 

173 


6. 797 

(. 024 ) 


.970 

(. 112 ) 


.950 
(. 088 ) 


6.743 

(. 027 ) 


.930 

(. 108 ) 


6.833 

(. 028 ) 


.940 
(. 700) 


M/G/l 
p=0. 90 


* 


0.99 

13.13C 


138,000 

69 


13.398 

(. 198 ) 


.675 

(. 243 ) 


.490 

(. 129 ) 


11.713 

(. 125 ) 


.820 

(. 288 ) 


12.174 

(. 119 ) 


.815 

(. 270 ) 




* 


0.99 

13.13C 


136,000 

17 


13. 240 

(. 202 ) 


.850 

(. 462 ) 


.725 

(. 274 ) 


11.467 

(. 125 ) 


.795 

(. 294 ) 


12.011 

(. 124 ) 


. 755 
(.277) 




* 


0.99 

13.13C 


128,000 

8 


13.011 

(. 200 ) 


.805 

(. 758 ) 


. 730 

(. 489 ) 


11.336 

(. 126 ) 


.755 

(. 286 ) 


11.762 

(. 115 ) 


.720 

(.200) 




t 


0. 999 
19.94 


1,386,000 

693 


20.693 

(. 324 ) 


. 700 

(. 161 ) 


.520 

(. 093 ) 


18.589 

(. 171 ) 


.890 

(. 211 ) 


19.063 

(. 157 ) 


.840 

(. 198 ) 




t 


0. 999 
19.94 


1,384,000 

173 


19.951 

(. 230 ) 


.890 

(. 294 ) 


. 800 

(. 170 ) 


18.188 

(. 159 ) 


.880 

(. 198 ) 


18.638 

(. 137 ) 


.880 

(. 185 ) 




t 


0.999 

19.94 


1,376,000 

86 


19.844 

(. 239 ) 


.820 
(. 385 ) 


. 780 

(. 238 ) 


18.053 

(. 176 ) 


.800 
(. 198 ) 


18.621 

(. 159 ) 


. 770 

(. 183 ) 




t 


0.999 

19.94 


5,536,000 

173 


20.152 

(. 146 ) 


.950 

(. 211 ) 


.920 

(. 156 ) 


19.511 

(. 138 ) 


.920 

(. 151 ) 


19.979 

(. 128 ) 


.950 

(. 142 ) 
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